Metal-insulator transition in the one- dimensional SU(N) Hubbard model 
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I. INTRODUCTION 

Although the metal-insulator transition has certainly 
been one of the most studied phenomenon in condensed- 
matter physics, it is only in recent years that impor- 
tant progress has been achieved. This is mainly due to 
careful experimental and numerical studies but also to 
the improvement of the theoretical tools [||-^). It has 
been proved extremely difficult to investigate the effect 
of strong correlations in dimensions greater than one, and 
it is only quite recently that, thanks to a new dynami- 
cal mean-field, our understanding has substantially pro- 
gressed M. For one-dimensional systems, the situation is 
rather different: There exist powerful analytical and nu- 
merical approaches at our disposal. Moreover, from the 
experimental point of view, the Mott-transition can be 
realized in organic conductors || and quantum wires Q . 
Therefore, one may expect to gain a lot of information 
on the physics of the metal-insulator transition. 

In one dimension, it has been recognized very rapidly 
that umklapp processes are at the heart of the problem. 
In the Abelian bosonization formalism, one can draw a 
general and consistent picture of the Mott-transition. In- 
deed, the charge properties are expected to be described, 
in the absence of umklapp contributions, by a Luttinger 
liquid with only two independent parameters: The charge 
velocity u c and the charge exponent K c that controls the 
decay of correlation functions. These quantities, which 
are non-universal, completely characterize the low-energy 
properties of a one-dimensional system . Within this 
framework, the effect of umklapp processes are investi- 
gated in perturbation theory, and one can write down 
an effective theory that describes the Mott-transition as 
well as a full description of the transport properties for 



any commensurate filling |s|]lOt . The only parameter that 
controls the transition is the (in general unknown) Lut- 
tinger charge exponent K c , and the transition is pre- 
dicted to be universal of the Kosterlitz-Thouless (KT) 
type. 

Most of the theoretical work in d = 1 focused on the 
properties of the standard SU(2) Hubbard model which, 
is known to be a Mott insulator at half-filling from its ex- 
act solution . Some extension of this model was con- 
sidered by introducing long-range hopping or finite-range 
interaction (nearest neighbor interaction for instance) [^J . 
In the present work, we study a most natural generaliza- 
tion of the usual Hubbard model: Instead of considering 
fcrmions with a two- valued spin index (with a SU(2) sym- 
metry) we generalize to the case of an arbitrary SU(N) 
spin index. Apart from the theoretical interest it is im- 
portant to emphasize that these additional degrees of 
freedom are realized physically through orbital degen- 
eracy as for example in Mn oxides ||. In this paper, 
we shall study the phase diagram of the one-dimensional 
SU(N) Hubbard model for repulsive interaction and at 
half filling corresponding to one "electron" per site. The 
Hamiltonian on a finite chain with L sites that we shall 
consider reads: 



n = -t 



L N 




(1) 



where the fermion annihilation operator of spin index 
a = 1,..,N at site i is denoted by Ci a and satisfies the 
canonical anticommutation relation: 



{c M ,c]j = <S Q iA 



(2) 
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The density of species a at the i site is denned by: 
n-ia = c la c ia- In the following, we shall consider that the 
nearest-neighbor hopping (t) and the on-site interaction 
(U) are positive. 

The Hamiltonian (]l|) is not exactly solvable by the 
Bethe ansatz for N > 2 and arbitrary U. It is however 
possible to solve the generalization of the Lieb-Wu Bethe 
ansatz equations for fermions carrying a SU(N) spin in- 
dex (HQ. The result is that for any N > 2, there 
exists a Mott-Hubbard transition from a metallic phase 
to an antiferromagnetic insulating phase at a finite value 
of the coupling U. The transition is found to be of first- 
order in contrast with the accepted view that the metal- 
insulator transition in one-dimensional systems should be 
of the KT type. The point is that a projection onto the 
subspace of states having at most two electrons at each 
site is crucial for the use of the Bethe ansatz approach. 
The other configurations are automatically excluded by 
the Pauli principle in the SU(2) Hubbard model whereas 
for N > 2 it is no longer the case. As a consequence, 
it is believed that the lattice model associated with the 
SU(N) generalization of Lieb-Wu Bethe ansatz equations 
should coincide with an integrable non-local version of the 
SU(N) Hubbard model @ @,P|- Although one natu- 
rally expects that the true SU(N) Hubbard model will 
share some properties with its non-local partner, in par- 
ticular the existence of a metallic phase at small enough 
U, the first-order character of the transition could take 
its origin in the nonlocality of the interaction. In any 
case, in order to study (|l|) one must abandon the exact 
Bethe ansatz approaches and resort to two powerful tech- 
niques available in one dimension: the bosonization and 
numerical approaches. As we shall show, none of these 
techniques is by itself sufficient to demonstrate the exis- 
tence of the Mott transition. Regarding bosonization, the 
mere existence of the metal-insulator transition -even in 
the simplest scenario of a KT phase transition- relies on 
the knowledge of U dependence of the Luttinger param- 
eter K c , a nonuniversal quantity which can only be com- 
puted in a perturbative expansion in U. In other words, 
bosonization cannot tell us whether a given lattice model 
will undergo a Mott-U transition. However, it defines a 
rich theoretical framework in which many qualitative and 
quantitative predictions are obtained. This provides an 
essential guide for the numerical investigation of a par- 
ticular lattice model. Regarding numerical investigations 
the situation is also not fully satisfactory. Beyond the ev- 
ident problem of memory and CPU time limitations, it 
is well-known that it is very difficult to characterize a 
KT phase transition. As we shall emphasize later, it is 
almost impossible to discriminate between the opening 
of a charge gap at U — and at a finite positive U, even 
when very accurate numerical data are at our disposal. 
The strategy employed in this work will consist in com- 
bining both approaches. Very strong evidences will be 
given in favor of a metal-insulator transition occuring at 
a finite positive value of the interaction U for N > 2. 



Various numerical methods can be used to study the 
ground-state properties of Hamiltonian (|l|). In exact di- 
agonalization methods the exact ground-state eigen- 
vector is calculated. Unfortunately, the rapid increase 
of the size of the Hilbert space restricts severely the at- 
tainable system sizes. In order to treat bigger systems 
two types of approach are at our disposal: The density 
matrix renormalization group (DMRG) method and the 
stochastic approaches. 

Since its discovery a few years ago the DMRG method 
has been extensively used for studying various one- 
dimensional systems and coupled chain problems (for a 
review, see Ref. J^L for a detailed presentation of the 
method, see Refs. DMRG is a very efficient real- 

space numerical renormalization-group (RG) approach. 
The fundamental point which makes the method suc- 
cessful is the way that "important" degrees of freedom 
are chosen at each RG iteration. Instead of keeping the 
lowest eigenstates of the RG block considered as isolated 
from the outside world (as it was usually done in pre- 
vious approaches), the states which are selected are the 
most probable eigenstates of the density matrix associ- 
ated with the block considered as a part of the whole 
system. The main error of DMRG is related to the finite 
number of states kept at each iteration of the algorithm. 
In order to get the exact property the extrapolation to an 
infinite number of states has to be performed. At least 
for ID and quasi-lD problems, and for systems having 
a small number of states per site, the errors obtained 
are small. Note also that DMRG works especially well 
when open boundary conditions are used. For periodic 
boundary conditions, errors are significantly larger. 

In this paper we use an alternative approach based on 
a stochastic sampling of the configuration space. Such 
approaches are referred to as quantum Monte Carlo 
(QMC) methods. There exists a large variety of QMC 
approaches. A first set of methods is defined within 
a finite-temperature framework (Path-Integral Monte 
Carlo, World-line Monte Carlo, etc- • •, see e.g. Ref. 
PH ). In these approaches, the main systematic error is 
the high-temperature approximation associated with the 
Trotter break-up (Trotter or short-time error). When in- 
terested in obtaining the zero-temperature properties the 
number of "time slices" to consider must be taken large 
and the computational effort becomes important. Prac- 
tical calculations have shown that the method is much 
less accurate than DMRG, at least for one-dimensional 
systems. In the second type of approaches used here, 
the stochastic sampling is directly defined within a zero- 
temperature framework. These methods are usually re- 
ferred to as Green's function Monte Carlo (GFMC) or 
projector Monte Carlo. For systems having a nodeless 
ground-state wave function as it is the case here, the 
GFMC method can be extremely powerful. The basic 
idea is to extract from a known trial wave function ipx its 
exact ground-state component ipQ. To do that an opera- 
tor G(Tt) acting as a filter is introduced. Statistical rules 
are defined in order to calculate stochastically the action 
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of the operator G on a given function. Apart from statis- 
tical fluctuations, the GFMC method is an exact method. 
It does not require an extrapolation to zero temperature 
as in finite-temperature schemes. In addition, there ex- 
ists a so-called zero- variance property for the energy: The 
better the trial wave function ipT is, the smaller the sta- 
tistical fluctuations are. In the limit of an exact wave 
function, the statistical fluctuations entirely disappear 
(zero- variance property). As an important consequence, 
by choosing a good enough trial wave function very accu- 
rate calculations can be performed (see, for example, Ref. 
fill). Note that, in contrast with DMRG, the efficiency 
of GFMC does not depend on the specific type of bound- 
ary conditions chosen and that the number of states per 
site is not a critical parameter of the simulation. Here, it 
is an important point since the SU(N) model displays 2 N 
states per site (for the SU(4) case treated here it gives 
16 states per site). 

In order to improve further the accuracy of the ap- 
proach we present a generalized version of the GFMC 
method in which the dynamics of the Monte Carlo pro- 
cess is partially integrated out. More precisely, we gener- 
alize an idea introduced by Trivedi and Ceperley in their 
GFMC study of the S=l/2 Heisenberg quantum antifer- 
romagnet jl9|. In the GFMC method the probability 
that the random walk remains a certain number of times 
in the same configuration is described by a Poisson dis- 
tribution. It is then possible to sample the correspond- 
ing "trapping time" from this distribution and to weight 
the expectations values according to it. As remarked by 
Trivedi and Ceperley, doing this can lead to a consider- 
able improvement in the simulation. This is particularly 
true when the wave function is localized (large U regime 
for our model, systems with deep potential wells, etc- • ■). 
Here, we show that the method can be improved fur- 
ther by integrating out exactly the time evolution associ- 
ated with this trapping phenomenon. Once this is done 
we are left with a random walk defined by an "escape 
transition probability" connecting non-identical configu- 
rations (the system never remains in the same configura- 
tion) and a modified branching term resulting from the 
time-integration. Note that introducing trapping times 
in averages helps a lot when optimizing the parameters 
of the trial wave function. Finally, we present an orig- 
inal method for computing the Luttinger liquid param- 
eters within a QMC scheme. We show that these pa- 
rameters can be obtained from a series of ground-state 
calculations of total energies of real -but not necessarily 
Hermitian- Hamiltonians. In this way we escape from the 
difficulty of calculating with QMC ground-state energies 
of the complex Hamiltonians resulting from the definition 
of the charge and spin stiffnesses. Although it is difficult 
to compare the efficiency of our generalized GFMC ap- 
proach with DMRG (since the quality of GFMC simula- 
tions is too much dependent on the quality of the trial 
wave function used) we believe that the accuracy of our 
results is comparable or even better to what can be done 
with DMRG. In any case, our data are sufficiently ac- 



curate to conclude to the existence of a metal-insulator 
phase transition in the model studied. 

Very recently, Beccaria et al. 20 have proposed a 
QMC algorithm based on the use of Poisson processes. 
Their approach contains similar ideas. However, in con- 
trast with the present approach no importance sampling 
is used and no integration of the Poisson dynamics is 
performed. It should also be noted that the use of Pois- 
son processes for describing the time evolution of systems 
trapped in some configuration is not restricted to quan- 
tum systems. Krauth and collaborators have proposed 
related ideas within the context of classical Monte Carlo 
simulations J2^Q. 

The organization of the paper is as follows. In section 
II, a bosonization approach of the SU(N) Hubbard model 
will be given. Some of the results has already been ob- 
tained by Affleck 0] whereas additionnal new ones will 
also be useful to compare with the numerical simulations. 
The purpose of section III is to give a presentation of the 
GFMC method together with our generalization based 
on the partial integration of the dynamics. The pratical 
implementations of the GFMC approach for the Hamilto- 
nian (|l|) will be discussed in section IV and the numerical 
results for N = 2, 3, 4 will be presented in section V. Fi- 
nally, section VI gives a summary of the work together 
with a comparison between the physical results obtained 
for the SU(N) Hubbard model and those corresponding 
to its nonlocal integrable version. In the Appendix we 
give some details of computation occuring in section II. 



II. THE SU(N) HUBBARD MODEL 

In this section, we shall use a bosonization approach 
(for recent reviews see Refs. piM) to study the SU(N) 
Hubbard model. Before doing that, let us first discuss 
the symmetries of the model. 

The Hamiltonian ([!]) has a U(l)(g> SU(N) symmetry: 



e l9 c 

UabC- 



ib 



(3) 



where the matrix U belongs to SU(N). These symmetries 
express the conservation of the charge and spin invariance 
under a SU(N) rotation. The associated generators are 
given by the following operators: 



with 



A/" = n ia 

i 

°i ~ c ia 2 ab c ib 



(4) 



(5) 



where the summation over repeated indexes (except for 
lattice indexes) is assumed in the following. In the latter 
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equation, the N 2 - 1 matrices T A are the generators of the 
Lie algebra of SU(N) in the fundamental representation. 
They satisfy the commutation relation: 



[T A ,T B ]=tf 



ABCq-C 



(0) 



f ABC being the structure constants of the Lie alge- 
bra and the generators are normalized according to 
Tr(T A T B ) = S AB /2. The conservation law associated 
with the U(l) symmetry allows to study the Hamilto- 
nian (Q) for a fixed density n. The Coulomb interaction 
can thus be rewritten, up to a constant, in terms of the 
SU(N) spin operator: 




UN 



A C A 



N+l 



S*S. 



where we have used the identity: 



rj-Arj-A 

1 ab 1 de 



bd 



1 

N" 



(7) 



(8) 



The relation (f7j) makes explicit the SU(N) invariance of 
the model. 

The Hamiltonian (|l|) is not exactly solvable by the 
Bethe ansatz for N>2 and arbitrary U, even if, as al- 
ready emphasized, some integrable non-local extension 
of (Q) with a SU(N) symmetry can be considered. The 
situation is simpler in the limit U — > oo and at half fill- 
ing (one "electron" per site), i.e. when kp — n/Nao 
(ao being the lattice spacing). In that case, it can be 
shown that (|l|) reduces to the SU(N) Heisenberg anti- 
ferromagnetic chain for which an exact solution is avail- 
able. As shown by Sutherland |25y , this latter model 
is critical with N — 1 massless bosonic modes with the 
same velocity. In the Conformal Field Theory (CFT) 
language, the central charge of the model in the infrared 
(IR) limit is c = N — 1 and using a non-Abelian bosoniza- 
tion of (Q) , Affleck [£3| identifies the nature of the critical 
theory in the spin sector as the SU(N)i Wess-Zumino- 
Novikov-Witten (WZNW) model. In the following, we 
shall present both non-Abelian and Abelian bosonization 
approaches of the SU(N) Hubbard model (0) at half fill- 
ing and give a number of results that will be essential for 
discussing the numerical data presented in section V. 



A. Continuum limit 

In the continuum limit, the spectrum around the two 
Fermi points ±kp is linearized and gives rise to left- 
moving fermions V'aL and right-moving fermions ip a R- 
In this low-energy procedure, the lattice fermion opera- 
tors Ci a are expressed in terms of these left-right moving 
fermions as: 

-% -> (X) ~ ^aR (X) e ik * X + lf> aL (x) e"**"", (9) 



where x = iaa. In this continuum limit, the non- 
interacting part of the Hamiltonian (]l|) corresponds to 
the Hamiltonian density of N free relativistic fermions: 



Ho 



-iv F {: ipl R d x ip aR : - : ipl L d x i/> a L :J (10) 



where the normal ordering :: with respect of the Fermi 
sea is assumed and the Fermi velocity vf is given by: 



vf = 2tan sin — . 

N 



(11) 



In the continuum limit, the SU(N) spin operator (^|) de- 
composes into an uniform and a 2k f contribution: 



a 



S A (a) ~ J A (x) + {c 2lkFX N A (x) + H.c.) (12) 



where the 2fcp contribution is given by: 

N A = ^ aL T A ^ bR (13) 
whereas the uniform part reads: J A — J R + J A with 

3R(L) = : ^aR{L)^bR(L) = • (14) 

These left-right SU(N) spin currents obey the following 
Operator Product Expansion (OPE) (see the Appendix): 



AD 



s 

Km Jr {L) (x) J b l) (y) ~ , , 



8ir 2 (x — yY 

^2^TT) J ^ ^ (15) 



fABC 



which shows that they satisfy the SU(N)i Kac-Moody 
(KM) algebra |Q,^6). In the same way, the total charge 
density J2 a nia rea ds in the continuum limit: 

J2 "*» «0 /2 (j° ( X ) + (^ 2lkFX ^aR (*) (x) + H.C. 

a 

(16) 

where J° = J% + J% and 

^R(L) =■ ^a R{ L)^aR(L) ■ (17) 

are the U(l) right and left charge currents. These cur- 
rents satisfy the OPE: 



hm J° (L) (x)J° (L) (y) 

x *y 



N 



4tt 2 (x — yY 



(18) 



and J belongs to the U(1)jv KM algebra. 

With these identifications, it is not difficult to show 
(see the Appendix) that the free part of the Hamiltonian 
( |iTj| ) can be expressed only in terms of spin and charge 
currents (the so-called Sugawara form): 



Ho — Hos + Hoc 



(19) 
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with 



and 



2itvf 
N + l 



<Jr j r. 



7TVf (■ tO-tO 



(20) 



(21) 



At the level of the free theory, spin and charge degrees 
of freedom decouple. The symmetry of the free Hamilto- 
nian 7io in the continuum limit is therefore enlarged to 
give: U(1) L ® SU(N) L ® U(l) fl (g) SU(N) jR , The Hamilto- 
nian Ho s is not hing bu t the Sugawara form of the SU(N)i 
WZNW model |£J,^6| . It is a conformaly invariant theory 
with central charge c = N — 1 (N — 1 massless bosons). 
The contribution 7i 0c describes the U(l) charge degrees 
of freedom and has central charge c = 1 (1 massless bo- 
son). 

The non-trivial part of the problem stems from the 
Coulomb interaction (Q). At sufficiently small U <<t, 
from Eq. (]?]), we see that its contribution will be given 
by the OPE: 

N 

V(x)=-Ua — lim S A (x + e) S A (x) . (22) 

V ' N+l e^O V / W \ I 

From Eq. @, there arc three contributions to V: 

V = V + V 2feF +V 4fcF . (23) 

The first term is the uniform k = component while 
the others contain oscillating factors e ±2lkpX and e ±4lkFX . 
Neglecting all oscillatory contributions, we are thus left 
with the uniform part Vo- Performing the necessary 
OPEs (see the Appendix), one finds that the total effec- 
tive low energy Hamiltonian density separates into two 
commuting charge and spin parts: 



H = Tic + T~i s 



(24) 



with 



n c = ^(: J&l : + : fifi :) + G c J°J° (25) 



and 



n s = 



27TW S 



j (: JUr ■■ + ■■ JlJ? ■■) + G s J£j£ (26) 



where the renormalized velocities are: 

Ua 



V s = V F 



2vr 



v c = v F + {N - 1) 



2tt 



(27) 



and the current-current couplings in the charge and the 
spin sectors are given by: 



G s = —2U etQ. 



Apart from a velocity renormalization, the effect of the 
Coulomb interaction is exhausted in the two marginal in- 
teractions in both charge and spin sectors. When U > 0, 
the spin current-current interaction is marginal irrele- 
vant. At the IR fixed point, G* = 0, the Hamiltonian 
in the spin sector is that of the SU(N)i WZNW model. 
On the other hand, the current-current interaction in 
the charge sector is exactly marginal since on can diago- 
nalize TL C with a Bogolioubov transformation to recover 
the Tomonaga-Luttinger Hamiltonian. Therefore, TL C de- 
scribes the line of fixed points of the Luttinger liquid. 

From the above analysis we conclude that the SU(N) 
Hubbard model at half filling is massless for small U > 0. 
The spin sector is described by the SU(N)i WZNW 
model while the charge sector is a Luttinger liquid with 
continuously varying exponents. The main point in the 
above analysis is the absence of umklapp terms which, 
when N —2, opens a gap in the charge sector for an 
infinitesimal value of the interaction. At this point it is 
worth recalling that the main approximation made in the 
above analysis is the omission of the oscillating contribu- 
tions V2k F and V4k F ■ This is a reasonable assumption as 
far as U is not too large. However one expects, on gen- 
eral grounds, that umklapp processes should contribute 
at sufficiently large U and that a Mott transition to an in- 
sulating phase should occur at a finite U c . Indeed, in the 
U — * oo limit, we have an insulating phase where the spin 
degrees of freedom are described by the SU(N) Heisen- 
berg antifcrromagnet. We shall return to this point later. 
For now let us concentrate on the properties of the metal- 
lic phase. 



B. The metallic phase 

At this point, we introduce N chiral bosonic fields 
l} aR(L)i a — (1, N), using the Abelian bosonization 



of Dirac fermions |2J] : 

^aR(L) 



exp (±iV47T(f> aR ( L) 



(29) 



where the bosonic fields satisfy the following commuta- 
tion relation [0 q ,r,0&l] = jfiab- The anticommutation 
between fermions with different spin indexes is realized 
through the presence of Klein factors (here Majorana 
fermions) n a with the following anticommutation rule: 
{K a ,Kb} — 25 a b- As in the N — 2 case, it is suitable 
to switch to a basis where the charge and spin degrees 
of freedom single out. To this end, let us introduce a 
charge bosonic field & c r(l) an d N — 1 spin bosonic fields 
^msR(L), m = (1, ...,N — 1) as follows: 



cR(L) 



+ - + M R ( L ) 



msR(L 



(28) 



) = 7=7 ==p= (0i + •■■ + <t>m - m<p m+1 ) R{L) . (30) 
y'm{m +1) 
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The transformation ( ^0| ) is canonical and preserves the 
bosonic commutation relations. The inverse transforma- 
tion is easily found to be: 



N-l 



IR(L) = -7=®cR(L) + 



IsR(L) 



5 afl(L) 



R(L) 



a- 1 



N 

N-l 



4> 



IsR(L) 



(a-l)sR(L) 



2,..., N-l 



4>NR(L) ~ —^^cR(L) 



N- 1 



N 



(31) 



In this new basis, the Hamiltonian density in the spin 
sector at the SU(N)i fixed point reads: 



N-l 



K = Y E ( : ( d ^ns) 2 : + : (d x Q ms f :) (32) 
where u s is the spin velocity at the fixed point and 

9ms = ^msi ~~ &msR- 



(33) 



This representation makes clear the fact that the central 
charge in the spin sector is indeed c = N — 1. 

Let us now concentrate on the charge sector. It is not 
difficult to show, using Eqs. (0, [|| and @ that the 
charge current expresses as: 



7° 



(34) 



Therefore, the Hamiltonian density ( p5|) in the charge 
sector reads: 



Vc 

2 



(d x $ c 



■■ (d x e c f 



+ (N-l)— d x $ cL d x $ cR 



(35) 



where we have introduced the total charge bosonic field: 
$ c = ®cR + tc£ and its dual: 6 C = <5> cL - $ c ii- The 



Hamiltonian (35) can be written in the Luttinger liquid 
form: 



(36) 



where the charge exponent K c and the renormalized 
charge velocity u c are given by: 



Kr 



1 



y/1 + (N - l)Uao/miF 

V F y/l + (N -l)Ua /TTV F . 



(37) 



The U dependence of the Luttinger parameters K c and 
u c given in the above expressions should not be taken 
too seriously. Indeed, the continuum limit approach is 
strictly speaking valid only provided U/t << 1. In this 
regime one has 



I- {N-l) 



Ua Q 



2irv F 

u c ~v F + {N-lf^. 



(38) 



The physically relevant question is now what happens 
for higher values of the interaction U. In the absence of 
umklapp terms, the accepted view is that the effect of 
interaction corresponds to a renormalization of the Lut- 
tinger parameters K c and u c as well as the spin velocity 
u s which have therefore to be thought as phenomenolog- 
ical parameters as the Landau coefficients in the Fermi 
liquid theory MM. These parameters completely char- 
acterize the low energy properties of the metallic phase 
as we shall see now. Let us first discuss the electronic 
Green's function defined by: 



G ab (x,T) = {iPI(x,t)MV,0)), 



(39) 



t being the imaginary time. This correlation function 
can be computed using Eqs. (|, fH and ©. After 
some calculations, one finds: 



G ab (x,T) 



2n 



1 



,/2 



X + U C T 

exp (ik F x) 



(iX + U C T) (IX + U S T) 

exp (—ikpx) 

+ 7 : , sl/JV / '■ ', -,1-l/JV 

( — IX + U C T) ( — IX + U S T) 

where the exponent a is given by: 
1 



2NK r 



{1-K C 



(40) 



(41) 



This allows to give an estimate of the single particle den- 
sity of states which is related to the electronic Green's 
function at x = 0: 



p(w) 



(42) 



Similarly, K c determines the singularity of the momen- 
tum distribution n a (k) around the Fermi point hp: 

n a (k) = n a {k F ) + Cte sgn (k - k F ) \k - k F \ a (43) 

and the momentum distribution function has a power law 
singularity at the Fermi level unlike a standard Fermi 
liquid. This anomalous power law behavior for any finite 
value of TV is inherent of a Luttinger liquid. 

The computation of the SU(N) spin-spin correlation 
function: 



G 



. AB 



(x,r) = (S A (x,t)S b (0,0)) 



(44) 



is more involved. It can be shown that the leading asymp- 
totics of this correlation function is given by the 2kp part: 



. AB 



(x,t) 



s:AB 



COS (2kpx) 



(x 2 + U 2 .T 2 ) Kc ^ N (x 2 + U 2 T 2 ) 



1-1/N ' 



(45) 



We deduce from the above correlation function the low 
temperature dependence of the NMR relaxation rate Ti : 



1 



j,2/N+2K c /N- 



(46) 



As seen, once the U dependence of the Luttinger pa- 
rameters u c , K c and u a is known, the low energy proper- 
ties of the metallic phase are entirely determined. These 
parameters are non-universal and cannot be obtained for 
arbitrary U by the continuum limit approach. Although 
K c < 1 when U > 0, one does not know its minimum 
value. It is only in the N — 2 case, that the Lut- 
tinger parameters can be extracted from the exact so- 
lution p7|-p9|]. When N >2 no exact solution is available 
and one has to use numerical computations to estimate 
these parameters. This will be done for the two cases 
N = 3 and N = 4 in section V. Before doing that, let 
us discuss about the Mott transtion that should occur in 
the problem for a finite critical value of the repulsion U 
for N > 2. 



C. The Mott transition 

The very difference between the N = 2 and N > 2 
cases lies in the fact that there is no umklapp term at half 
filling in the bare Hamiltonian in the continuum limit. 
The reason for this is that these terms came with oscillat- 
ing factors and were omitted for small value of the repul- 
sion. However, in the RG strategy one has to look at the 
stability of the Luttinger fixed line and any operator that 
is compatible with the symmetry of the problem should 
be taken into account: they will be generated during the 
renormalization process. In our problem, the important 
symmetries are the SU(N) spin rotation invariance, chiral 
invariance and translation invariance. 

From Eqs. M, |29|) and (||(J), one easily finds that under 
a translation by one lattice site, the charge field <I> C is 
shifted according to: 



4\ 



<S>r 



(47) 



Therefore one can add any operator in the charge sector 
that is invariant under the transformation (^7j) and will 
be necessary generated by higher order in perturbation 
theory. The operator with the smallest scaling dimension 
that is invariant under (47) is: 



n 



umklapp 



-G u : cos 



(48) 



Other operators, with higher scaling dimensions, that 
couple spin and charge degrees of freedom may also be 
included. This is the reason why one cannot exclude the 
possibility of a charge density wave (CDW) instability. 
For instance, such processes are present in the extended 
SU(2) Hubbard model at half filling Q. Alhough it re- 
quires some formal proof, we expect that, due to the fact 
that in the present model the interaction is local in the 
density, the leading umklapp contribution should only 
affect the charge sector. We have checked that this is 
indeed true for the particular cases, N — 3 and N = 4 
pl[ . We have shown indeed by perturbation theory that 
the oscillating contributions V2k F and V^ F generate 6kp 
and Akp processes for TV = 3 and N = 4 respectively. Up 
to irrelevant operators, the only contribution we found 
is precisely ( |48|) with N = 3 and N = 4. In any case 
in what follows, we shall thus make the hypothesis, first 
made by Affleck [p3f , that all the effects of high energy 
processes are exhausted by ( fl8| ) for the general SU(N) 
case. Consequently, the effective Hamiltonian density in 
the spin sector is still given by the SU(N)i WZNW model 
and the effective Hamiltonian in the charge sector is now: 



K c : (d x @ c y 



-G u : cos ( V4ttN$ 



and C 



(49) 



Rescaling the fields as $ c 

the Hamiltonian ( [49] ) in the charge sector becomes the 
Hamiltonian of the sine-Gordon model: 



H c = 



(: {d x $> c f : + : (d x Q c ) 
—G u : cos 



(50) 



Since the scaling dimension of the cosine term in ( |50| ) is 
A u = NK C , we deduce that a gap opens in the charge 
sector when: 



2 

N' 



(51) 



On the other hand, when K c < 2/N, the umklapp term 
is irrelevant and the system remains in the metallic phase 
described in the preceeding subsection. Therefore, as U 
increases, K c will decrease from 1 at U — to K c = 2/N 
at a critical value of the interaction U c where a Mott tran- 
sition to an insulating phase occurs. Within this scheme, 
the phase transition is expected to be of the KT type. 
Of course when U > U c , K c vanishes so that the jump is 
1 — 2/N and is universal. The present approach cannot 
give an accurate value of U c . However, one can get a 
rough estimate of U c using Eqs. (O, ^) and (|5l|): 



Uc 
t 



7T 7V2 

2 



4 . 7T 

sm — . 

N - 1 N 



(52) 
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In the insulating phase, the charge field <I> C is locked 
in a special well of the sine-Gordon model ( [sof ) and the 
leading asymptotics of the SU(N) spin-spin correlation 
functions is now: 



A ab (x,t)^X 1 6 



AB cos(2kpx) 

(x 2 + U 2 r 2)l-l/iV 



(53) 



where Ai is a non-universal constant stemming from the 
charge degrees of freedom. One recovers the result previ- 
ously derived by Affleck (23). The NMR relaxation rate 
behaves now as l/(TiT ) - T 2 /^ 2 . Finally, let us note 
that there are other harmonics 4kF,6kp.. in the SU(N) 
spin density (jl^) that will be generated by higher or- 
ders in perturbation theory. Together with the uniform 
contribution with scaling dimension 1, these terms will 
give subleading power law contributions in the SU(N) 
spin-spin correlation function ( p3[ ) . These operators cor- 
respond to the primary fields of SU(N)i WZNW trans- 
forming to other representation of SU(N) than the fun- 
damental one. One should recall that for the SU(N)i 
WZNW, there are N — 1 primary fields |2(J . A primary 
field O (a = 1, .., N — 1) of SU(N)i transforms according 
to the & th basic representation of SU(N) (Young tableau 
with a boxes and a single column) and has scaling dimen- 
sion: A a = a(N — a)/N. We thus expect the following 
asymptotics for A AB with some non-universal constants 
(Aa): 



A AB (x,r) 



KAB 



8tt 2 



f 



f 



{u s t — ix) {u s t + ix) 



N-l 



+5 AB V A g cos ( 2afe ^) (54) 
a tt (x 2 +u 2 rr- a/N 

up to logarithmic contributions originating from the 
marginal irrelevant current-current interaction in the 
spin sector p2[ . 

We end this subsection by giving the low-temperature 
expression of the uniform susceptibility x an( i the spe- 
cific heat of the SU(N) Hubbard model in the insulating 
antiferromagnetic phase. The continuum density that de- 
scribes the behavior of the SU(N) spins degrees of free- 
dom in a uniform magnetic field H is given by: 

iV-l 

Uh = t E ( : ( d ^ns) 2 ■■ + ■■ (d x e ms y 



N-l 



-H J™ (55) 



m— 1 

where we have neglected the marginally irrelevant 
current-current interaction. In Eq. (pq), we have con- 
sidered a uniform magnetic field along the diagonal 
T m ,(m = l,..,N - 1) generators of SU(N) that span 
the Cartan subalgebra of SU(N). According to our nor- 
malization convention, they can be written in N x N 
diagonal matrices as follows: 



^J2m (to + 1) 



diag(l,l,...,-m,0,..,0) (56) 



with to = 1, .., N — I and —to is located on the m + 1 
element of the diagonal. Using the bosonization corre- 
spondence ( ^9| ) and the canonical transformation (|3~l|), 
the total density Hamiltonian ( |55| ) in a magnetic field 
can be written as: 

N-l 

Hh = Y ^ ( S * $ ™) 2 : + : ( d - e " ls f :) 



H 



N-l 



/2vr 



Doing the substitution: 

d $ — > d $ 



H 



2iru s 



(57) 



(58) 



we obtain the expression of the uniform susceptibility of 
the SU(N) Heisenberg antiferromagnet: 



X 



N-l 



(59) 



which is nothing but N — 1 times the uniform suscep- 
tibility of the SU(2) Heisenberg antiferromagnet. This 
result is easy to understand since the critical theory in 
the spin sector corresponds to N — 1 decoupled mass- 
less bosonic modes. Finally, using the general formula 
of the specific heat at low temperatures for a conformaly 
invariant theory [p3|, one has for the SU(N) Heisenberg 
antiferromagnet : 



ir (N-l) 
3m s 



T. 



(60) 



Before closing this section, it is important to empha- 
size that the Mott transition expected in the bosoniza- 
tion approach relies on the full expression of K C (U) as 
function of the interaction. However, one should stress 
that this parameter cannot be obtained for arbitrary U 
within this approach and only in the weak coupling limit 
U <t where the model is in its metallic phase. To con- 
clude in favour of the existence of a Mott transition for 
a finite value of the Coulomb interaction, one has thus 
to compute K C (U) of the lattice model by an indepen- 
dent approach. Since the SU(N) Hubbard model with 
N > 2 is not exactly soluble, one cannot determine the 
expression K C (U) by the Bethe ansatz as for the stan- 
dard Hubbard model p7|~p9| . We shall thus compute the 
value K C (U) of the lattice model using very acurate nu- 
merical calculations based on QMC methods described 
in the next section. In section V, we shall then com- 
pare the numerical results with the predictions of the 
bosonization approach given in this section to conclude 
on the existence of a Mott transition in the model. 
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III. THE NUMERICAL APPROACH 

In this section we present our improved zero- 
temperature Green's function Monte Carlo method used 
for computing ground-state properties. In the first part 
a sketchy but self-contained presentation of the basic 
GFMC method is given. In addition to introducing our 
notations for the next part, this section will enable the 
interested reader to understand all the practical aspects 
of the method. The second part is devoted to the pre- 
sentation of the generalized GFMC method itself. 



A. Green's function Monte Carlo 

As already noticed in the introduction the basic idea 
of the GFMC method is to extract from a known trial 
wave function \ipT > the exact ground-state component 
\tpo >. To do that an operator G(H) acting as a fil- 
ter is introduced. For continuum problems standard 
choices are G(H) = exp(— tTL) (Diffusion Monte Carlo) 
or G{H) = 1/(1 + t(E - H)) (Green's function Monte 
Carlo). For a lattice problem or any model with a fi- 
nite number of states (finite matrix) a natural choice to 
consider is 



G(H) = 1-t(H- Et) 



(61) 



where r plays the role of a time-step (a positive constant) 
and Et is some reference energy. If r is chosen sufficiently 
small and | ipx > has a non-zero overlap with the ground- 
state, the exact ground-state is filtered out as follows 



lim G(H) \tpT >~ |-0o > 



(62) 



This result is easily obtained by expanding \ipx > within 
the complete set of eigenstates of 7i. 

In Monte Carlo schemes, successive applications of 
the operator G(H) on |?/>t > are done using proba- 
bilistic rules. These rules are implemented in configu- 
ration space where the trial wave function and matrix 
elements of 7i are easily evaluated. In what follows we 
shall denote by | i) an arbitrary configuration of the 
system. To give an example, in actual calculations pre- 
sented below we consider | i) = \i^> > ■ ■ ■ \i( N > > with 
| >= \nia, • • • , riL a > where L is the number of sites, 
a the SU(N) color index, and rii a the occupation number 
of site i (riia — or 1) for the species a. 

In this work Hamiltonians considered are of the form: 



U = T+ V 



(63) 



where T is the kinetic term (a non-diagonal operator) 
and V is a (diagonal) potential term. For fermions in 
one dimension it is known that by choosing a suitable 
labelling of the sites, matrix elements of the kinetic term 
can all be made negative 



<i\T\j> <0 (i^j). 



(64) 



A most important consequence of this property is that 
the exact ground-state has a constant sign. In other 
words, simulations presented here are free of the sign 
problem. 

Let us now introduce the following transition proba- 
bility 



P^(r)=V>T(jXi| [l-r(H-E L )]\i) 



1 



(65) 



where ipT(i) are the components of the vector \tpT >, 
V>t(z) =< i\i>T >, and El is a diagonal operator called 
the local energy and defined as follows 



< i\E L \j >= S i:j E L (i) 



with 



El{i) 



< j\n\ip T > 

< i\ip T > 



(66) 



Note the important relation associated with the defini- 
tion of the local energy: 



(H - E L ) | Vt> = 0. 



(67) 



To define a transition probability Pj_>j must fulfill the 
two following conditions. First, the sum over final states, 
P%^j, must be equal to one. Here, this is true as a di- 
rect consequence of (|67|). Second, -Pi— >j must be positive. 
To see this and for later use, let us distinguish between 
the cases i = j and i =/= j. 
For i = j we have 



Pi^i(r) =1 + tTz(») 



(68) 



where T L (i) = E L {i) - H vl . Using (|63j), T L (i) can be 
rewritten as 



<i\T\<p T > 
< i\ipr > '' 



(69) 



Tl{i) is called the local kinetic energy. Because of (f34|) it 
is a negative quantity and the transition probability can 
be made positive by taking r sufficiently small. More 
precisely, the time-step must verify 



< r < Mmi[-l/T L (i 



(70) 



Note that the upper bound is a non-zero quantity for a 
finite system. 

On the other hand, when i ^ j we have 



(71) 



a positive expression since ipT(i) is chosen to be positive 
and off-diagonal terms Hij are negative (Eq. (|64)). 

Using expressions ( |68| ) and ( [71] ) for the transition prob- 
ability random walks in configuration space can be gen- 
erated. By averaging over configurations, statistical es- 
timates for various quantities can be defined. A first 
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important example is the calculation of the variational 
energy associated with \ipT > (variational Monte Carlo). 
The variational energy is defined as 



E v (ifr) = 



Here, it is rewritten as 



< 1pT\T-C\i>T > 



1 - 

E v (ip T ) = lim -V E L {i) =« E L » 

K— >oo A A — ' 



(p) 



(72) 



(73) 



where << ... >>(p) is the stochastic average over con- 
figurations \i > generated using the transition probabil- 
ity P, K being the number of configurations calculated. 
Equation ( [73] ) holds because ?/>t(z) 2 is the stationary den- 
sity of the stochastic process, that is 



5> t « 2 ^(t) = Mj) 2 Vj. 



(74) 



This property is directly verified by using expressions 
©and®. 

As already pointed out, the estimate of the ex- 
act energy is based on the stochastic calculation of 
[1 - t{H - E T )] n \ipT >, Eq. @. Introducing between 
each operator in the product the decomposition of the 
identity over the basis set, 1 = J2i N >< *l) anc ^ making 
use of the definition of the transition probability, (65), 
we get the following path integral representation 



p-i 



[1 - t(H - E T )f | ipr) = ^( J °) 2 II P > 



fc=0 



p-i 

n 

k=0 



"ikik+l 



ip) 



where the weights tu !3 are defined as follows 
(* | [1 - t(H - E T )] | j) 



W; 



[i\[l-r(H-E L )}\j) 



(75) 



(76) 



or more explicitly, 



From 



Wij = 1 i ^ j, 

_ 1 - t{H u - Et) . _ . 
hl l-T(H u -E L (i)) l - J - 

the exact energy can be obtained as 



E 



lim 



< ifr\H[l - t{H - Et)] \ih > 



(77) 



(78) 



p^oo < il> T \[l - T (H - E T )] \iI>t > 
which is rewritten here in terms of stochastic averages as 



E 



lim 

P^oo 



p-i P-i 
<< E L (i P ) Y[ w ihik+1 »(p) / << II »(P) 



k=0 



k=0 



(79) 



In order to compute the averages appearing in that ex- 
pression two strategies can be employed. First, formula 
([79|) can be directly used as it stands: Paths are gener- 
ated using the transition probability and the local energy 
at each step is weighted by the quantity W — J| fc Wi k i k+1 . 
This approach where the number of configurations is kept 
fixed and the weights are carried out along trajectories is 
usually referred to as the Pure Diffusion or Green's func- 
tion Monte Carlo method. For extended systems such as 
those considered here, this approach is not optimal. In- 
deed, it is important to sample less frequently regions of 
configuration space where the total weight is small and 
to accumulate statistics where it is large. To realize this, 
a birth-death process (or branching process) associated 
with the local weight is introduced. In practice, it con- 
sists in adding to the standard stochastic move defined 
by the transition probability, a new step in which the 
current configuration is destroyed or copied a number of 
times proportional to the local weight. Denoting m%j the 
number of copies (multiplicity) of the state j, we take 



rriij = m.t(wij + rf) 



(80) 



where int(a;) denotes the integer part of x, and 77 a uni- 
form random number on (0, 1). Adding a branching pro- 
cess can be viewed as sampling with a generalized tran- 
sition probability P£_>j (r) defined as 



P*^.(r) = Pi^(r) Wij 
M3){j I [l-r(H-E T )]\i) 



(81) 



Of course, the normalization is not constant (the popula- 
tion fluctuates) and P* is not a genuine transition prob- 
ability. However, we can still define a stationary density 
for it. From Eq. ( |si|) we see that the stationary condi- 
tion is obtained when Et is chosen to be the exact energy 
Eq, and that the density is ipr(i)ipo{i)- Accordingly, by 
using a stabilized population of configurations the exact 
energy may be now obtained as 



En = « Et >> 



(P,w) ■ 



(82) 



Note the use of an additional subscript w in the average 
to recall the presence of the branching process. 

At this point, we shall not expand further the method. 
For more details regarding the implementation of GFMC 
to lattice systems the interested reader is referred to Refs. 
p4] , p5 19 3^]. Let us just emphasize on two important as- 
pects. First, there exists a so-called zero-variance prop- 
erty for the energy: The better the trial wave function 
ipT is, the smaller the statistical fluctuations are. In the 
limit of an exact wave function for which the local en- 
ergy is a constant, fluctuations entirely disappear (zero 
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variance). From this important remark follows that in 
any QMC method, it is crucial to optimize as much as 
possible the trial wave function used. Of course, in prac- 
tice, a compromise between the complexity of the wave 
function and the gain in reduction of variance has to be 
found. 

Once a good trial wave function has been chosen, the 
only room left for improvement is the implementation 
of the dynamical process itself. In the algorithm pre- 
sented here the only dynamical parameter which can be 
adjusted is the time-step r. In a configuration \i > asso- 
ciated with a small local kinetic energy Tl(i), the system 
remains in this configuration a relatively large time and 
a large value of t is necessary to help the system to es- 
cape from it. Unfortunately, because of the constraint 
( pq ) (P,-_»j must be positive) configurations with a high 
local kinetic energy impose a small value of r. In order 
to circumvent this difficulty, we propose to integrate out 
exactly the time evolution of the system when trapped 
in a given configuration. This idea is developped in the 
next section. 



B. GFMC and Poisson processes 

Consider the probability that the system remains in a 
given configuration i a number of times equal to n. It is 
given by 

Vi(n) = P(it =i,T;...;i n = i,T; i n +i ^ i, r) = 

[P^(r)ni-P^(r)] (83) 

Vi(n) defines a normalized discrete Poisson distribution. 
In terms of the local kinetic energy it can be rewritten as 



1-1 



Vi(n) = -rT L {i) exp [nln(l + tT l {%))\ 



(84) 



where the integer n runs from zero to infinity. To describe 
transitions towards states j different from i we introduce 
the following escape transition probability 



Pi, 



l-Pi^(r) 



3 ^ i- 



(85) 



Using Eqs. ( p8p and fl69|) Pi~>j is rewritten in the most 
explicit form 



P, 



n lj ip T {i) 



(86) 



Note that this transition probability is positive, normal- 
ized, and independent of the time-step r. Now, by using 
both probabilities, Vi(n) and Pi^j, the path integral rep- 
resentation of G(H.) P \4>t >, formula ([75]), can be rewrit- 
ten as 

[1 - t(H - E T )] P | tpr) = J2 ^( l °) 2 

(i,n)eC P 



[\\P lk (n k )P k ^ k+l \P tl (ni)\{ 



1 



k=0 



k=0 



H) (87) 



where the sum is performed over the set of all families 
of states (iq . . . ii) with multiplicities (uq . . . n{) verifying 

J2k=o( n k+l)+ n l = P- In a given family successive states 
are different and the variable rik represents the number 
of times the system remains in configuration i^. The set 
of all families is denoted Cp and an arbitrary element 
is written (i,n) = (i$ . . . ii,no . . .ni). Since off-diagonal 
weights are equal to one, Eq. ([77]), a shortened notation 
for the diagonal weights, Wi = wu, has been introduced. 

Now, let us remark that the time-step r plays a role 
in the path integral formula (^) only when the system 
is trapped into a given configuration. Indeed, both the 
escape probability P and the off-diagonal weight Wij are 
independent of r. As an important consequence the limit 
t — > and P — > oo with Pt — t can be done exactly. In 
this limit the discrete Poisson process Vi(n) defined in 
( gjj ) converges to a continuous Poissonian distribution 
for the variable 9 = nr 



1 



9/6 



In this formula 9i represents the average time spent in 
configuration i. In what follows we shall refer to it as the 
average trapping time, its expression is 



Si = -l/T L (i). 



(89) 



The fact that 9i is inversely proportional to the local ki- 
netic energy is explained as follows. When the kinetic 
energy is small the system is almost blocked in its config- 
uration and 9 is large. In contrast, when a large kinetic 
energy is available, the system can escape easily from its 
current configuration and 9 is small. As already remarked 
the escape transition probability is independent of t and 
is therefore not affected by the zero-time-step limit. Fi- 
nally, after exponentiating the product of weights, the 
path integral can be rewritten in the form 



-t(H-E T ) 



1pl 



E 



dS t 



o ■ • 



l-l 



[n^(^)P fc ^ +1 ]^(^)e-SU ^(^fe)^' 



•) 



1 



fe=0 



(90) 



with the constraint that the trapping times verify 
E fc= o °k = t. 

In order to compute ground-state properties the limit 
t — > oo must be performed, Eq. (p2[). In this limit the 

constraint ^2 k=0 Ok = t can be relaxed and, quite re- 
markably, integrations over the Poisson distributions for 
the different trapping times can be performed. For large 
enough time t we obtain 
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l-l 



e ^(*o) 2 n^-u- +1 n 



-i i 



H) (91) 



k—0 k—Q 

where the new integrated weights w are found to be 

T L [i) 



ET — 



(92) 



In the same way as before the exact energy can be ob- 
tained as 

Eq = hm — — . (93) 

In terms of stochastic averages it gives 
Eq = lim 

I — >oo 



<< E L (ii)Q ll JT w lk » {P) i « e k Yl w lk » 



fc=0 



fe=0 



(P) 



(94) 



where configurations are generated using the escape tran- 
sition probability P. 

As in the standard approach it is preferable to simulate 
the weights via a branching process. Here also, the ref- 
erence energy Et stabilizing the population is given by 
the exact energy, Eq. The new stationary density writes 



7r(i) ~ ipT{i)i>o(i)/0i 



(95) 



up to an immaterial normalization constant. Finally, our 
estimator for Eq is 



Eq = 



« 9iE L (i) » ( p^) 



<< Vi » 



(96) 



(P,w) 



where configurations are generated using P and branched 
with w. Note that the variational energy can be recovered 
by removing the branching process (w = 1) 



« 6iE L (i) »,p 



« Vi » { p ) 



(97) 



A. Hard-core boson Hamiltonian 

The Hamiltonian considered here is the one- 
dimensional SU(N) Hubbard model described by (Q). 
Simulations are performed for a finite ring of length L. 
In one dimension the sites can be labelled in such way 
that the hopping term connects only sites represented by 
consecutive integers. As a consequence no fermion sign 
appears, except eventually when a fermion crosses the 
boundary (l^IorL->l). By choosing either periodic 
or antiperiodic boundary conditions this sign can always 
be absorbed and our model ([!]) becomes equivalent to a 
model made up with hard-core bosons and described by 



L N 



U = t J2 E c m^a + H.c. + V - E(E ™») 2 (98) 



i— 1 a— 1 



where cf a creates a hard-core boson of color a on site i, 
Uia is the occupation number rii a = cf a Ci a , and c^ +ia = 



B. Trial wave function 

As already emphasized a most important aspect of any 
Monte Carlo scheme is the choice of a good trial wave 
function. To guide our choice, let us consider the exact 
solution at U = 0. In this case the ground-state is ob- 
tained by filling N independent Fermi seas consisting of 
planes waves with momenta k n = 2Tm/L (n = 0, ±1, . . .). 
For a given type of fermion, the ground-state can be writ- 
ten as a Vandermonde determinant J37| and the following 
expression for the ground-state is obtained 



<=°(ii, 
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n 

Kl' 



sin — (it 



ii>)] 



(99) 



where i±, . . . ,ip are the positions of the P fermions on 
the chain, ik = 1, . . . , L. In terms of occupation numbers 
the solution can be rewritten as 



t nAQn 



<j>(ni,. . .,n L ) 



(100) 



where the matrix Aq of size (L x L) is given by 



IV. COMPUTATIONAL DETAILS 

In this section some important aspects of the practical 
implementation of the GFMC approach to the SU(N) 
Hubbard model are presented. 



Ao(i, i') = In | sin[— (i ~ i') 



(101) 



Note that (100) and (101) describe a system of particles 
interacting via a logarithmic potential (ID Log-gas). The 
exact ground-state wavefunction of the complete SU(N) 
model at U = is simply obt ained by writting the prod- 
uct of the N wavefunctions (100) associated with each 
color. 

When the Coulomb interaction is switched on, we have 
chosen to take the same functional form as before for ipT 



12 



t nA\jn 



(102) 



Here, Ajj is an arbitrary matrix of size (NL x NL). Tak- 
ing into account the translational and SU(N) symmetries, 
at most L + 2 independent variational parameters can be 
defined. In all GFMC calculations presented in this paper 
the entire set of parameters has been systematically op- 
timized. To do that, we have generalized the correlated 
sampling method of Umrigar et al. |38| along the lines 
presented in the preceding section. To be more precise, 
the set of configurations used to calculate the quantities 
to be minimized (variational energy or variance of H, 
see Ref. Q) are generated using the escape transition 
probability and weighted with the corresponding average 
trapping times. Doing this, the effective number of con- 
figurations is increased and the optimization process is 
facilitated. We have found that that large numbers of 
parameters can be easily optimized. 



C. O(L) algorithm 

In the occupation-number representation the numeri- 
cal effort for calculating the trial wave function iprin) is 
of order 0(L 2 ). To evaluate the local energy the Hamil- 
tonian has to be applied to the vector \ipT >■ Since a 
given configuration | n) is connected by H to about O(L) 
states, the total computational cost per Monte Carlo step 
is about 0(L 3 ). In fact, this cost can be reduced to 0(L). 
To do that, we introduce the following set of 2NL + 1 
variables 



f - \-f- a - nA v n 
(n, nu,n ) = [n, A v n, — - — ). 



(103) 



Using this representation, the wave function is given by 
e n ° . Configurations connected by the Hamiltonian differ 
from each other by removing a particle of a given color a 
on a site i and putting it on a neighboring site j. In the 
occupation-number language it corresponds to add one to 
the component ja and remove one to the component ia 
of vector n. For convenience let us introduce the vector 
<j( M ) whose components are zero except the component 
ia which is equal to one. Using the new variables just 
defined we have 

(n,n^,n ) -> 



(n + 5^ - P a \nu + Au8 {ia) - A v 5 {ia \ 



n 



t0(ia) _ $J°))A u (P a ) - <fa'«)) 



nu 



8 Ua) )) 
(104) 



In the simulation the set of new variables is stored for 
each configuration. At each Monte Carlo step they are 
reactualized using (104). Finally, the numerical effort is 
limited to O(L). 



V. RESULTS 

Let us now present the results for the SU(2)-, SU(3)- 
and SU(4) Hubbard models. SU(2) results have been 
obtained by solving numerically the Lieb-Wu equations 
p| . Other results have been obtained with the GFMC 
method presented in the previous section. In all calcula- 
tions we have set t = 1. 



A. Charge gaps 

The finite-size charge gap A c (N e , L) is defined as 

E {N e + 1,L) + E {N e — 1,L) — 2E (N e ,L) (105) 

where Eq(N b , L) is the total ground-state energy of a ring 
of length L with N e electrons. In this expression N e ± 1 
means that a fermion of an arbitrary color is added to 
or removed from the system. Denoting N the number 
of colors, calculations are done for a number of fcrmions 
of each color equal to L/N, and therefore for a total 
density n = N e /L equal to one. In order to get the exact 
charge gap the limit L — > oo must be performed. As usual 
this is done by calculating charge gaps for different sizes 
and extrapolating to infinity. Here, SU(3) and SU(4) 
calculations have been done for L = 9, 12, 18, 27 and L = 
8, 16, 24, 32, respectively. The finite-size gaps have been 
found to converge almost linearly as a function of the 
inverse of the size. Accordingly, the limit L — > oo of 
the gap has been obtained from a fit of the data with a 
linear or quadratic function of 1/L. Figure [l] presents the 
charge gaps obtained for N = 2, 3, 4 as a function of the 
Coulomb interaction U. 
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Fig. 1. Charge gaps as a function of the interaction U for 
the SU(2), SU(3), and SU(4) Hubbard models. The values of 
the gaps have been extrapolated to L — * oo (see text). 
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A first important remark concerns the quality of the 
Monte Carlo simulations. As it can be seen in Fig. [j], 
the error bars on the different gaps are quite small. A 
typical value is about 0.001. Errors are small because 
total energies are calculated with a very high level of ac- 
curacy. For example, for the SU(4)-modcl with L = 32 
and U = 0.5, we get £ (32,32) = -52.13056(15) for a 
total number of elementary Monte Carlo steps equal to 
8.10 7 . Clearly, the relative error of about 3.10" 6 is very 
small. In the large [/-regime where the trial wave func- 
tion is not expected to be as good as for small U, we 
still get excellent results. For example, for U = 4.5 we 
get £ (32,32) = -23.7118(13) (1.610 s MC steps) with 
a relative error of about 6.10 -5 . Using the standard 
GFMC method (presented in Sec. III. A) we get, for 
U = 0.5, E Q {32, 32) = -52.13050(40) and, for U = 4.5, 
-E (32,32) = -23.7210(110) (in both cases the maximum 
time-step allowed has been chosen, see Eq. ([70|)). The 
improvment resulting from the new approach, particu- 
larly at large U, is noticeable. Finally, using the ap- 
proach of Trivedi and Ceperley Jl9| (introduction of the 
Poisson process but no integration in time) we get for 
U = 0.5 £ (32,32) = -52.13041(22) and for U = 4.5, 
E {32, 32) = -23.7121(30). These results illustrate the 
improvment resulting from the time integration. 

Having at our disposal such accurate results we can 
try to find out whether or not a gap opens for a non-zero 
value of U. Considering only continuous transitions, two 
scenarios are possible. A first possibility is to open a gap 
for any non-zero value of U. In that case we write the 
gap versus U as follows 



A c = Cexp(-G/U). 



(106) 



A second scenario consists in looking for the existence of 
a KT-type transition at a finite value U c of the Coulomb 
interaction. In that case the gap is written as 



A c = C K t exp (- 



G 



KT 



VU-u c 



(107) 



for U > U c and zero otherwise. The thre e set s of results 
have been fitted either using Eqs. (106) or (|107[). The 



fitting procedure used is a standard one, based on the 
minimization of a chi-square type function including sta- 
tistical errors. Our most important conclusion is that all 
sets of data can be correctly represented within our small 
stati stica l errors either using the gapful representation, 
Eq. (106), or using a KT scenario, Eq. (107) with a not 
too large value of U c . For example, using Eq. (106) pos- 
sible representations are (C = 25.313, G = 11.318), (C = 
274.634, G = 26.745), and (C = 515.649, G = 32.755), 
for N — 2,3, and 4, respectively. Although no clear 
physical content can be given to the magnitude of coeffi- 
cients, it is nevertheless s atisf actory to verify that in the 
case of SU(2), the gapful (106) leads to not too large val- 
ues for the coefficients. This should be contrasted with 
the SU(3) and SU(4) cases for which the parameters are 
important. Within a KT scenario all data can also be 



very well fitted. In the case of SU(2) where we know 
for sure that no KT transition exists, the 'critical value' 
issued from our fits ranges from to about 0.5. For 
example, a possible representation is given by (Ckt — 
541.310, G KT = 11.053, and U c = 0.384). For the Sub- 
model accurate representations can be obtained with a 
value of U c ranging from to about 2.3 For U c = 2.2 
(the value we shall propose later for the critical value) 
we get: (C KT = 45.050, G KT = 6.567, and U c = 2.2). 
For SU(4) the interval is larger. Allowed values range 
from to about 2.9. For U c — 2.8 (our proposed value, 
see below) we get: (C K t = 17.889, G K t = 5.144, and 
U c = 2.8). In contrast with the gapful representation, it 
should be noted that coefficients are now much larger for 
the SU(2) model than for the SU(3) and SU(4) models. 

In conclusion, using accurate values of the gaps no con- 
clusions can be reasonably drawn about the existence or 
not of a KT-type transition at a finite value of U. Nu- 
merical evidences based on other quantities are therefore 
called for (see next sections). From the fitting of our 
data the only conclusion we are allowed to draw is that 
a KT-transition is only possible within the range (0,2.3) 
for SU(3) and within the range (0,2.9) for SU(4). In ad- 
dition to this, if such a transition actually occurs in both 
models, we should expect a difference for the critical val- 
ues given by: U C [SU{4)] - U C [SU(3)} ~ 0.5-0.6 (see Fig. 



B. Spin gap 

The spin gap is defined as the change in ground-state 
energy produced when destroying a fermion of a given 
color and creating a fermion of a different color (in the 
SU(2) case it consists in flipping one spin). Note that in 
this process the charge number is kept fixed. For a finite 
system we have 

A s (A^ e , L) = E (N e + l,L)- E (N e -1,L) (108) 

where N e ± 1 involves an arbitrary pair of electrons of 
different colors. 

For the SU(2) case the system is known from the exact 
solution to be gapless for an arbitrary value of the inter- 
action strength U. For a number of colors greater than 2, 
it is an open question. This is an important point since 
the existence of a gapful regime would very likely indicate 
the existence of a coupling between spin and charge de- 
grees of freedom. In all calculations performed for N=3 
and 4, and for a coupling constant U ranging from very 
small to very large values (up to U — 10) no evidence 
for the existence of such a gap has been found. Thus, 
it can be quite safely concluded that the spin sector of 
SU(N) N = 2, 3,4 is gapless for an arbitrary interaction 
in full agreement with the bosonization prediction. To 
illustrate this point we present in Fig. || a typical behav- 
ior for the spin gap of SU(3) as a function of 1/L at the 
relatively large value U = 4.5 (at least two times greater 
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than the maximal value expected for U c in the charge 
sector). The behavior of the gap is essentially linear and 
extrapolation to the origin leads to a vanishing gap. 




Fig. 2. Spin gap as a function of l/L for the SU(3) Hubbard 
model at U = 4.5. The solid line is a linear fit of the data. 



C. Luttinger liquid parameters 

In this section we present calculations of the Luttinger 
liquid parameters u c and K c . For that we shall make use 
of their relations with the compressibility k and charge 
stiffness D c of the system. For a model with N colors 
(SU(N)) we have the following relations 



N 



and 



TTU C 2 

KTl = — 

Kr 2 



D r = Nu r K r 



(109) 



(110) 



where n — N e /L (N e total number of electrons) is the 
electron density. The compressibility k is defined as the 
second derivative of the ground-state energy Eq with re- 
spect to the density of particles 



1 d 2 E 



L dn 2 



(111) 



A convenient finite-size approximation of the compress- 
ibility is 

L 

tv — 77 

N 2 

E (N e + N,L) + E (N e — N,L) — 2E (N e ,L) 
N 2 

(112) 

where N e ± N in Eq means that N fermions -one of each 
color- are added to or removed from the system. 
The charge stiffness is given by 



D c = 



7T d 2 E 



L dip 2 



(113) 



<p=0 



where ip is a charge twist in the system. This charge twist 
is imposed by introducing the following twisted boundary 
conditions 



= f C 



(114) 



for an arbitrary site i and color a. 

By calculating with GFMC total ground-st ate e nergies 
for different numbers of electrons, formula (112) allows 
a direct calculation of the compressibility. In contrast, 
the GFMC calculation of the charge stiffness is more 
tricky due to the presence of a complex hopping term 
at the boundary. To circumvent this difficulty we resort 
to the second-order perturbation-theory expression of the 
charge stiffness. We have 



„ 7T 

D c = - 
L 



< —T > -2 



< fc|J|0 > 
Ek — E Q 



(115) 



12( C i+la C ia 



H.c.) is the kinetic-energy 
i+ia c ia - H.c.) is the paramag- 



where T 

operator, J = — it~Y^{c, 
netic current operator, < • • • > denoting the expectation 
value in the ground-state, all quantities being evaluated 
at ip = 0. To evaluate the kinetic term we make use of 
the Hellman-Feynman theorem: < T >= Eq ~ U^f-. In 
practice, the following finite-difference expression is used 



< T >= E — U 



E (U + 6U) - E(U - 
26U 



6U)\ 



(116) 



with 5U small enough to make higher-order contributions 
negligible. 

The second-order part of formula ( |115| ) can be re- 
interpreted back as the second-derivative of the total 
ground-state energy of a new Hamiltonian consisting of 
the original Hamiltonian plus a perturbation associated 
with the flux operator J. This leads to the following 
relation 



E 

k=£0 



< fcjJjO > 

Eq — Ek 



1 d 2 E (\) 

2 OX 2 



(117) 



where Eq is the ground-state energy of the new Hamilto- 
nian defined by 

H = -(t + A) Y^(ct +laCia ) - (t - A) ^ktla^a) + V(U) 

ia ia 

(118) 

and V(U) is the po tent ial part of the problem. Using 
formulas (117) and (118) the charge stiffness can now be 
obtained from a series of GFMC ground-state calcula- 
tions of total energies of reai Hamiltonians (more pre- 
cisely, E , Eq(5U), and E (-SU) for H, and E (X) for 
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H, Eq. flllSD ). It should be emphasized that the new 
Hamiltonian H is real but not symmetric: Left-moving 
and right-moving electrons do not have the same velocity. 
Of course, such a property is easily implemented within 
a QMC framework. 

Figures ||§|]^|||],H present the Luttinger parameters 
u c and K c for the SU(2), SU(3), and SU(4) Hubbard 
models as a function of the interaction U and for dif- 
ferent sizes L. For the SU(2) model, parameters have 
been obtained by computing ground-state energies issued 
from the standard Lieb-Wu equ ations (computation of 
the compressibility, formula (112)) and from their gener- 
alization to the case of twisted boundary conditions as 
presented by Shastry and Sutherland |39| (computation 
of the charge stiffness, formula ( |113| )). For the SU(3) and 
SU(4) models we have followed the general route just pre- 
sented above. 

A first striking result emerging from the figures is the 
strong qualitative differences between the general behav- 
ior of Luttinger parameters of the SU(2) model on the one 
hand, and of the SU(3) and SU(4) models, on the other 
hand. Let us first have a look at the charge velocity u c . 




Fig. 3. u c as a function of U for the SU(2) Hubbard model. 
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as a function of U for the SU(4) Hubbard model. 



Fig. 4. u c as a function of U for the SU(3) Hubbard model. 



In the SU(2) case the charge velocity has been cal- 
culated for various values of U and for the sizes L — 
6, 10, 14, 18, and 22. Results are presented in Fig. ||. 
The upper curve corresponds to L = 6, the lower one to 
the maximum size, L = 22. In between, the curves are 
ordered according to the magnitude of L. For a given size 
L, the charge velocity is found to decrease as a function 
of U. For a given U, u c also decreases as a function of 
the size L. Such a behavior is quite typical of a gapped 
system in which collective charge excitations are damped 
away. In the limit of large sizes, the charge velocity is ex- 
pected to vanish for a non-zero value of the interaction. 

The charge velocities of the SU(3) model, Fig. [|, and 
of the SU(4) model, Fig. [| display a very similar be- 
havior which is dramatically different from the one ob- 
served for SU(2). Starting from their free value at U = 
(u c = \/3 and u c = \[2 for SU(3) and SU(4), respec- 
tively), they increase as a function of U with a finite 
slope at the origin. After some critical value of U both 
velocities go down quite rapidly. In the first part of the 
curves (small and intermediate values of U) the charge 
velocity is found to converge quite rapidly as a function 
of the size. All curves presented cannot be distinguished 
within statistical errors. Although the calculations pre- 
sented here are limited to systems with a maximum size 
of L = 27 (SU(3)) or L = 32 (SU(4)) some preliminary 
calculations at larger sizes strongly suggest that the val- 
ues plotted are indeed converged. Such results strongly 
support the existence of a gapless phase for the SU(3) 
and SU(4) models. At larger values of U the situation 
is rather different. The charge velocities decrease quite 
rapidly both as a function of U and as a function of L. 
This behavior indicates the existence of a gapped phase. 
In order to be more quantitative let us have a look at 
the value of the slope at the origin. The theoretical pre- 
diction can be obtained from Eqs. (j38[). For SU(3) cal- 
culations have been done for the sizes L=9,18, and 32. 
The slope at the origin is found to be 0.32(1), 0.32(1), and 
0.33(2) for L=9,18, and 27, respectively. Theses results 
are in perfect agreement with the theoretical prediction 
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of — ~ 0.318. For the SU(4) model calculations have 
been done for the sizes L=16,24, and 32. The slope at 
the origin is found to be 0.46(1),0.47(1), and 0.45(2) for 
L=16,24, and 32, respectively. Here also, the results are 
in perfect agreement with the theoretical prediction of 
~ 0.477. Let us now consider our results for K c . 
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Fig. 6. K c as a function of U for the SU(2) Hubbard model. 




Fig. 7. K c as a function of U for the SU(3) Hubbard model. 




Here also, there exists a common behavior for the cases 
SU(3) and SU(4) and a different one for SU(2). In the 
latter case, Fig. ra, K c decreases either as a function of 
U or as a function of the size. The slope at the origin, 
U = 0, is essentially zero and K c is expected to vanish 
at large sizes, except, of course, in the free case. Once 
again, this behavior is typical of a gapped system. In 
the two other cases, the situation is rather different. In 
the same way as for the charge velocity, two regimes can 
be distinguished, see Figs. and ||. At small and inter- 
mediate U, the values of K c are found to be very well 
converged within statistical errors as a function of the 
size L. The curve is smooth with a finite slope at the 
origin. In the second regime corresponding to larger val- 
ues of U the curves K c versus U go down as a function 
of the size. Clearly, this latter regime corresponds to a 
gapped phase. Having nearly exact values of K c up to 
some critical value U c for SU(3) and SU(4), the next log- 
ical step consists in comparing these values to the predic- 
tions of bosonization. A first important prediction was 
the opening of a gap in the charge sector for a value of 
K c equal to 2/N, Eq. (|5l|). In Fig. [?] corresponding to 
the SU(3) dashed line has been drawn at the value 

K c = 2/3. The intersection of this line with the curves of 
K c appears at about U c ~ 2.2. A most remarkable result 
is that this value of U is both consistent with the critical 
value extracted from the calculation of the charge gaps, 
Fig. [ij but also with the fact that it lies in the domain 
of U where the values of K c begin not to converge as a 
function of the size (a fact usually interpreted as resulting 
from the existence of a finite correlation length) . A very 
similar situation is obtained in the SU(4) case. Using the 
same type of arguments, U c is found to be around 2.8. 
When studying charge gaps we had observed a difference 
of U c , Fig. |, between SU(3) and SU(4) of between 0.5 
and 0.6. This is in very good agreement with what is 
found here from independent data on K c . A second pre- 
diction which can be tested is the estimate of the value 
of U c itself. Formula ([52]) gives 
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Fig. 8. K c as a function of U for the SU(4) Hubbard model. 



For N = 3 and N = 4 one gets U c = 3.40 and U c = 4.44, 
respectively. As already pointed out, these estimates 
must be considered with caution. However, it should 
give the correct trend as a function of N. Here, if we 
look at the ratio U C [SU(4)]/U C [SU(3)] we get about 1.31 
from the theoretical estimate and about 1.27 from our 
data. The agreement is excellent. Another point which 
can be checked is the value of the slope at the origin. 
For the SU(3) case, it is found to be -0.18(1),-0.19(1), 
and -0.19(2) for L=9,18, and 27, respectively. Theses 
results are in very good agreement with the theoretical 
prediction of — —0.183 given by Eq. (|38|). For 

SU(4) we find a slope of -0.31(1),-0.33(1), and -0.32(2) 
for L=16,24, and 32, respectively. Theses results are 
also in total agreement with the theoretical prediction 
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of 



-0.337. 



Finally, it can be very useful for interested readers to 
give some compact and accurate representations of the 
Luttinger parameters K c and function of U. For 

both parameters a minimal representation we may think 
of (see section II. B) is 



1 



(119) 



u c = Vpy/l + UiU + u 2 U 2 . 
For SU(3) we obtain 

fci = 0.33452, k 2 = 0.08789 



ux = 0.37929, u 2 = -0.025509. 

Note that these values are not too far from the bare values 
corresponding to Eqs. @, fc? = u\ = — ~ 0.36755, 



uQ 

"-2 



= 0. 



For SU(4) we obtain 

fei = 0.62065, fc 2 = 0.12298 

U! = 0.71486, u 2 = -0.052705 

to compare to the bare values given by k® = u\ = — — = 
0.675237, k% = ul = 0. 

As already discussed we have found no evidences for 
the opening of a spin gap in the case of the SU(3) and 
SU(4) models. In other words, the system remains crit- 
ical with respect to the spin degrees of freedom for any 
value of the interaction. For these models the slope at 
the origin is predicted to be equal to — ^ ~ —0.159 (Eq. 
©). Once again, this value has been recovered using 
our numerical data. To compute the spin velocity we 
have used the formula expressing the spin gap as a func- 
tion of the size for a critical system |28| 



A,(jV e ,£) 
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(120) 



For SU(3) and SU(4) we get for the slope -0.18(2) and 
-0.18(3), respectively, in very good agreement with the 
theoretical prediction. 

A final piece of information which can be extracted 
from our data is related to the way the total ground- 
state energy converges to its asymptotic value. To be 
more precise, it is known that the ground-state energy 
per site, eo(L), of a Luttinger liquid is expected to behave 
as follows [E8| 



e (L) ~ e (+oo) - — ^ 



(121) 



where ■ ui denotes the total velocity associated with all 
critical excitations. In the free case, N degres of freedom 



are critical, and the total velocity is equal to Nvp. When 
the interaction is turned on, it is possible to follow the 
evolution of the total velocity as a function of U. This has 
been done for the SU(3) model. Taking our data for the 
sizes L=9, 18, and 27 the ground-s tate energy has been 
fitted with a form adapted to Eq. (121), eo = a — b/L 2 . 
From this fit an effective number of critical modes can be 
defined as follows 



Neff = 



ttvf 



The result is presented in Fig 




Fig. 9. Effective number of critical modes as a function of U 
for the SU(3) Hubbard model. 



Although the transition is not as sharp as for the Lut- 
tinger parameters, the loss of one critical mode (passing 
from 3 to 2) is clearly seen when U varies from zero to 
infinity. A similar curve may be obtained for the SU(4) 
case. 



VI. CONCLUDING REMARKS 

In this work, we have studied the SU(N) generalization 
of the one-dimensional Hubbard model for repulsive in- 
teraction at half-filling. Using a combination of bosoniza- 
tion and QMC results, we have clearly shown that the 
SU(N) Hubbard model for N > 2 behaves very differently 
from the SU(2) case. Strong numerical and theoretical 
evidences have been given in favor of a Mott transition, 
between a metallic and an insulating phase, occuring for 
a finite value of the Coulomb repulsion U c > for N > 2. 

The picture emerging from the bosonization approach 
consists in a spin-charge separation at low energy. The 
spin degrees of freedom are critical for arbitrary U and 
described by the SU(N)i WZNW model with a central 
charge c = N — 1 (N — 1 gapless bosonic modes) . The 
effective theory associated with the charge degrees of 
freedom corresponds to a sine-Gordon model at (3 2 — 
4ttNK c (U). For a small value of the Coulomb interac- 
tion U, the interaction is irrelevant. The charge sector 
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is then critical and described by a massless bosonic field. 
In this weak coupling phase, the system is metallic with 
anomalous power law behaviors in the physical quan- 
tities typical of a Luttinger liquid. For a finite value 
of the interaction U c such that K C (U C ) = 2/N, a KT 
phase transition to an insulating phase is expected in the 
bosonization approach. In this strong-coupling phase, 
the charge bosonic field becomes locked and the infinite 
discrete symmetry related to the periodicity of the 
potential of the sine-Gordon model is spontaneously bro- 
ken. The only degrees of freedom that remain critical 
in this strong coupling phase are the N — 1 spin modes 
and after integrating out the massive charge degrees of 
freedom, the low-energy theory of the model corresponds 
to the SU(N) Heisenberg antiferromagnet. 

Very accurate numerical simulations based on a gener- 
alization of the GFMC method and fully optimized trial 
wave functions have been performed to obtain the spin 
and charge gaps, and the Luttinger liquid parameters 
as a function of the Coulomb interaction for the SU(2), 
SU(3), and SU(4) Hubbard models. A metal-insulator 
phase transition at a finite value U c is clearly seen for 
SU(3) (U c ~ 2.2) and SU(4) (U c ~ 2.8) in contrast with 
the standard SU(2) case. In addition all the results ob- 
tained for N — 3 and N — 4 are fully consistent with 
the theoretical framework drawn in section II. This pro- 
vides an accurate test of the bosonization approach to 
the SU(N) Hubbard model for small and large values of 
U. It is therefore natural to expect that the physical pic- 
ture emerging from the two cases studied here can be ex- 
tended to arbitrary values of N. Thus one may conclude 
that the occurrence at a finite value of the interaction of 
a Mott transition of the KT type is generic in the SU(N) 
Hubbard model for N > 2 at half-filling. In addition, it 
should be emphasized that the calculations of the Lut- 
tinger parameters K c and u c presented in Sec. II. B are 
of very good quality (in particular they are converged as 
a function of the size) and thus provide an accurate char- 
acterization of the low-energy properties of the metallic 
phase of the SU(3) and SU(4) Hubbard models. 

Let us now compare our results with the exact solution 
of the integrable model based on the SU(N) generaliza- 
tion of the Lieb-Wu Bethe ansatz equations . As dis- 
cussed in the introduction, an exact solution of an SU(N) 
generalization of the Hubbard model is available. Al- 
though the underlying lattice Hamiltonian of the model is 
not known, it involves very likely long-range interactions 
that dynamically exclude three-electron configurations. 
The question that naturally arises is whether the physics 
described by the latter model is similar, when N> 2, to 
that of the lattice SU(N) Hubbard model that we have 
studied in this paper. At half-filling, the SU(N) inte- 
grable model undergoes a first-order phase-transition, 
as one varies U , from a metallic to an insulating phase 
p3[ . This is in disagreement with the KT transition pre- 
dicted by our analysis. In the metallic phase the inte- 
grable model is a Luttinger liquid for every N [ ^3|j4l| 
with the same physical properties as those obtained by 



the bosonization approach for the SU(N) Hubbard model. 
However, the charge stiffness K c obtained from the Bethe 
ansatz equations varies between 1/N and 1 as U de- 
creases from U c to Jl|^i|]. The value at the transi- 
tion (K c = 1 /N) is thus two times larger than the value 
obtained for the SU(N) Hubbard model. This clearly 
confirms that the integrable model differs from the lat- 
tice SU(N) Hubbard model in the charge sector. As al- 
ready pointed out, this difference should result from the 
presence of non-local interactions in the lattice model as- 
sociated with the integrable SU(N) model. 

Regarding perspectives, it is clearly of interest to fur- 
ther explore the phase diagram of the SU(N) Hubbard 
model: case of an attractive interaction, dependence 
on the filling, etc... For an attractive interaction at 
half-filling, bosonization predicts that a phase transition 
should also occur as \U\ varies. For incommensurate fill- 
ings, it is easy to see, within the bosonization framework, 
that the system is a Luttinger liquid for arbitrary TV and 
positive U where the leading asymptotics of the electronic 
Green's function and spin-spin correlation coincide with 
those computed in the metallic phase. The situation is 
less clear for commensurate fillings Uf = nn/ (Nao) (N/n 
being an integer). In the bosonization approach, a gap 
opens in the charge sector for K c = 2n 2 /N. The ex- 
istence of a Mott transition for commensurate fillings 
clearly requires the full knowledge of K c (U,n) of the 
lattice model. Some preliminary calculations show that 
there is a very special commensurate filling, n — N/2, 
where no Mott transition exists and for which the charge 
and spin degrees of freedom are massive for N > 2 and 
arbitrary U 40 1. 

Let us end by noting a very interesting connection 
bewteen the metal-insulator transition predicted in the 
SU(N) Hubbard model and the existence of plateaux in 
magnetization curves of spin ladders under a strong mag- 
netic field [^2|-^4| . Using the Jordan- Wigner transforma- 
tion, one can indeed interpret the SU(N) Hubbard model 
as a N-leg S=l/2 XY spin ladder in an uniform magnetic 
field along the z-axis and coupled in a symmetric way by 
Ising interaction. The relation between the Fermi mo- 
menta and the magnetization < M > (normalized such 
that the saturation value is ±1) is: kp = < M > 

)/(2ao). The Mott transition found in this work for the 
SU(N) Hubbard model at half-filling corresponds to the 
appearance of plateaux at < M >= (N — 2)/N in the 
magnetization curves of the previous N-leg XY spin lad- 
der. Moreover, the existence of a Mott transition for the 
SU(N) Hubbard model at commensurate filling will give 
additional plateaux located at < M >= (N — 2n)/N in 
the magnetization curves of the corresponding spin lad- 
der. 
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APPENDIX: 

In this Appendix, we give some details of computations 
to establish the separation of spin and charge (|5|) at the 
Hamiltonian level in the continuum limit of the SU(N) 
Hubbard model and fix the expressions of u c s and G CjS 
given by Eqs. @ M. 



1. Sugawara form of the free Hamiltonian 

To begin with, we shall recall some basic things on the 
SU(N) non-Abelian bosonization (for a review see Refs. 
[Hll|6)). As seen in section IIA, the chiral SU(N) spin 
current 3^ L can be expressed in terms of N left-right 
moving fermions ifj a R,L- 



(Al) 



The left (respectively right)-moving fermions are holo- 
morphic (respectively anti-holomorphic) fields of the 
complex coordinate (z — r + ix, r being the imaginary 
time): ijj a L(z), ipaR(z)- These fields are defined by the 
following OPEs: 

V"L ( z ) ^bL M 



2ir (z — uj) 

+ ■ ipl L ihL ■ M + (z-uj): dipl L ip bL : (uj) + .. 



^aR ( Z ) ^bR (W) 



Sab 



2tt (z — uj) 



+ '■ WaR^bR ■ M + (z-U)) : dipl R ip bR : (uj) + 



(A2) 



with d — d u , d — do and there are no singularities in the 
OPE when one does the fusion of two operators belonging 
to different sectors. 

Let us now consider the OPE between two left SU(N) 
spin currents for instance: 

Jl CO Jl M =: ^aL^bL : (z) : ^jJ^eL : M 

= T^Al (*) ^l («) (z) ^ dL M . (A3) 

Using the OPEs ([A|), the commutation relation (|^), and 
the normalization of the generators of the SU(N) Lie al- 
gebra, one obtains: 



Jl (z) JE (u) 



+ - 



7T 2 (z — Ul) 
: fABC 



2tt (z — ui) 

In the same way, we find for the right spin current 

6 AB 



(A4) 



Jg(z)Jg(uj) 



+ 



IT 2 (z — Uj) 
: fABC 

JR{-). 



2tt (z — uj) 



(A5) 



Evaluating these OPE at equal time, one recovers the 
OPE (15) showing that Jr L are SU(N) 1 spin current. 
With the same procedure, one can compute the OPE 
between the charge current J R L using its definition (|lj 
in terms of the underlying fermions: 



Jl(z)Jl(uj) 
J r(z)J rW 



N 



47T 2 (Z - Uj) 

N 



47T 2 (Z 



(A6) 



so that the charge current Jr L belongs to the U(1)jv 
KM algebra. 

The next step is to obtain the Sugawara form (^(j, [2l]) 
of the free part of the Hamiltonian (Hq). Let us consider, 
for instance, the left sector of the theory since we shall 
obtain the same result for the right part with the substi- 
tution: L — ► R, (z,w) — ► (z,uj) and d — ► 8 . We need 
now the following OPE for the spin sector: 

j£ (z) Jt (uj) =: i>t L T a ii> bL : (z) : 4 L T d A e ^ eL : (uj) 

= ~ (s ae S bd - ^5 ab 5 de j i/J aL (z) tjj eL (uj) i> bL (z) ip dL (uj) (A7) 



where we have used the relation (||). Using QA2| ) and 
keeping also the first regular terms in the fusion, we get: 



8tt 2 (z - ujY 



N ■ 



2N 



N 2 



2-kN 

Therefore, one obtains 
AM 



- : ^ f aL^aL^bL^l L ■ (uj) 



(A8) 



ST A rrA 

• Jl Jl •— 



2N 



- ■ ll>t L 1paLlpbLll>i L 
~ ■ ^aL d ^aL ■ 



2ttN 



(A9) 



In the same way, we obtain for the left charge current: 
: JlJl := - : ^aL^bL^L : ™ ■ ihfyaL ■■ • (A10) 

IT 

One can eliminate the four fermions terms by considering 
the following combination: 



7T 

N 



JlJl 



2tt 



N + 1 



Jr A Ji 



:^ aL d^ aL :. (All) 



Since one ha s dipa L — —id x ip a L within our convention, 
the identity ( All ), the so-called Sugawara form, states 
that the free Hamiltonian of N relativistic left-moving 
fermions can be written only as a function of left current- 
current terms. In the right part, we have also a similar 
identity: 
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N ' JrJr 



2tt 

'n + i 



st A rrA 

Jr. Jr 



-i : 4>i R d x iPaR : ■ (A12) 



Collecting all terms, we finalyobtain the Sugawara form 
of the free Hamiltonian 7io ( |l0|) : 



1 ( : ^ aR dx^aR ■ - ■ 1pl L d x 1pal 



j-f {■ JrJr 



-Jill ■■) 



2lT 



N- 



r (: JrJi 



A a- A 
R 



j£j£ 



). (A13) 



2. Sugawara form of the SU(N) Hubbard 
Hamiltonian 

We shall now investigate the effect of the Hubbard in- 
teraction in the continuum limit to fix the expressions 
( p7| ) and (p8j) of the velocities (u C;S ) and the coupling 
constants (Gc jS ). Using the continuum description of the 
SU(N) spin density (12|), the interacting part (0) is given 
by dropping all oscil 



12|) , the interacting part 
atory contributions: 



Vn 



Ua N 
' N + 1 



J A J J 



: Af A ::N A ^ : + 
: N M :: N A :) 



(A14) 



The OPE between the 2kp p arts of the spin density can 
be computed using (13) and (A2) as in the previous sub- 
section. We find up to constant terms: 



{z, z) : M A ^ : (w, Q) + : : {z, z) : J\f A : (u, w) 
N 2 - 1 z - uj 
2ttN z-uj 
N 2 - 1 z - Q 



WaL^aL ■ (w) 



: ^i R dipaR ■ (w) 



2?r7V z- 

•4>t L tpaLtPb R ^bR ■ (W,0>) 



A/ 



1pl L tpbLlpl R llJaR ■ . 



Using Eqs. (A9, AlC) and similar equations in the right 
sector together with the definition of the charge current 
(|l7|), we end with: 



J A J A +N A N A ^ +N A ^N A ~ - 



iV 2 



= 

and 

W s = 



/7O <t0 . 

JrJr ■ 



JlJl 



2ttv s 
N+l 



ST A ST A 

Jr. Jr. 



rrA rrA 

Jl Jl 



Gc 3131 (A18) 

f~l rrA rrA 

(j s Jr Jl ■ 

(A19) 



The renormalized velocities are given by: 
Ua 



vf 



2tt 

(A-l)^ 

Z7T 



(A20) 



whereas the current-current couplings in the charge and 
the spin sectors write: 



G r = — — — U a 



G s 



N 
-2Ua . 



(A21) 



(A15) [7; 



2N 2 



(: J°r.J°r 



1 

N 1 
N 2 



ST A ST A 

Jr. Jr 



j£j£ 



,N+1 
' N 



JkJl 



-70 /70 

N 2 JrJl- 



(Al^ 



As a consequence, the continuum limit of the SU(N) Hub- 
bard model at half filling exhibits the spin-charge sepa- 
ration: 



TL = Tic + Ti s 



(A17) 
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